#include "ScreeningCalculos.h"


void ScreeningCalculos::calcular()
{
    #pragma omp single
    {
        if (modelo->SPRList.Begin!=NULL)valorDaVez=modelo->SPRList.Begin->foiAcessado;
        CriaMatriz();
        for (int j = 1; j <= Cols; j++)
        {
            for (int i = 1; i <= Rows; i++)
            {
                //Pego cada elemento de uma célula
                for (SPHERE* s = Cell[i][j]; s != NULL; s = s->CellNext)
                {
                    //E testo com todos os demais da própria célula e das células vizinhas
                    for (SPHERE* t= s->CellNext;    t != NULL; t = t->CellNext) Forca->ForcaResSS(s,t);
                    for (SPHERE* t= Cell[i-1][j-1]; t != NULL; t = t->CellNext) Forca->ForcaResSS(s,t);
                    for (SPHERE* t= Cell[i-1][j  ]; t != NULL; t = t->CellNext) Forca->ForcaResSS(s,t);
                    for (SPHERE* t= Cell[i-1][j+1]; t != NULL; t = t->CellNext) Forca->ForcaResSS(s,t);
                    for (SPHERE* t= Cell[i  ][j-1]; t != NULL; t = t->CellNext) Forca->ForcaResSS(s,t);
                    for (LINE* l = Lines->Begin; l != NULL; l = l->Next) Forca->ForcaResSL( s, l );
                    if (rigido)
                    {
                        if (s->R[0] < BoxMinH + s->Rad && s->V[0] < 0.0) { s->R[0] = BoxMinH + s->Rad; s->V[0] = -s->V[0]; }
                        if (s->R[1] < BoxMinV + s->Rad && s->V[1] < 0.0) { s->R[1] = BoxMinV + s->Rad; s->V[1] = -s->V[1]; }
                        if (s->R[0] > BoxMaxH - s->Rad && s->V[0] > 0.0) { s->R[0] = BoxMaxH - s->Rad; s->V[0] = -s->V[0]; }
                        if (s->R[1] > BoxMaxV- s->Rad && s->V[1] > 0.0) { s->R[1] = BoxMaxV - s->Rad; s->V[1] = -s->V[1]; }
                    }
                    else
                    {
                        if (s->R[0] < BoxMinH + s->Rad
                        || s->R[1] < BoxMinV + s->Rad
                        || s->R[0] > BoxMaxH - s->Rad
                        || s->R[1] > BoxMaxV - s->Rad)
                        {
                            s = modelo->Delete(s);
                            continue;
                        }
                    }
                    s->Force( s->Mass, gravidade );
                    if (MaterialOceanico)
                    {
                        REAL invDens;
                        REAL raioQuad;
                        REAL normaV;
                        REAL mult;
                        invDens = -1.0 * DensidadeMatOc;
                        raioQuad = s->Rad * s->Rad;
                        normaV = NORM(s->V);
                        mult = invDens*raioQuad;
                        mult *= normaV;
                        s->Force( mult, s->V );// Equação do arrasto (Lord-Rayleigh)
                    }
                    checkSprings(s);
                }
            }
        }
    }
}

void ScreeningCalculos::checkSprings(SPHERE *s)
{
    for (SPRING* k = s->Springs; k != NULL; )
    {
        SPRING* prox=k->SNext(s);
        //Evitar computar a mesma mola quando encontrar a outra esfera conectada
        if ((k->foiAcessado)==valorDaVez)
        {
            //printf("real ");
            REAL Dist, Dist2, F;
            REAL R[2], V[2];

            k->foiAcessado=!valorDaVez;
            k->Translate();

            R[0] = k->R2[0] - k->R1[0];
            R[1] = k->R2[1] - k->R1[1];
            Dist2 = DOT( R, R );
            Dist = sqrt(Dist2);

            V[0] = k->S2->V[0] - k->S1->V[0];
            V[1] = k->S2->V[1] - k->S1->V[1];
            F = ( k->K * (k->L0 * Dist - Dist2) - k->C * DOT( V, R ) ) / Dist2;

            k->S1->Force( -F, R );
            k->S2->Force( F, R );

            k->S1->Torque( -F * (PERP_DOT( R, k->R1 ) - PERP_DOT( R, k->S1->R )) );
            k->S2->Torque( F * (PERP_DOT( R, k->R2 ) - PERP_DOT( R, k->S2->R )) );

            k->MeasureStress( Dist );
            if (k->Stress > k->Mat->Sigma_E) k->L0 += k->L0 * (k->Stress - k->Mat->Sigma_E) / k->Mat->Young;
            else if (-k->Stress > k->Mat->Sigma_E) k->L0 += k->L0 * (k->Stress + k->Mat->Sigma_E) / k->Mat->Young;

            if (Dist > k->FractureT || Dist < k->FractureC)
            {
                k=modelo->Delete(k);
            }
        }
        k=prox;
    }
}

